Self-organized criticality model for brain plasticity 
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Networks of living neurons exhibit an avalanche mode of activity, experimentally found in organ- 
otypic cultures. Here we present a model based on self-organized criticality and taking into account 
brain plasticity, which is able to reproduce the spectrum of electroencephalograms (EEG). The model 
consists in an electrical network with threshold firing and activity-dependent synapse strenghts. The 
system exhibits an avalanche activity power law distributed. The analysis of the power spectra of 
the electrical signal reproduces very robustly the power law behaviour with the exponent 0.8, ex- 
perimentally measured in EEG spectra. The same value of the exponent is found on small-world 
lattices and for leaky neurons, indicating that universality holds for a wide class of brain models. 

PACS numbers: 05.65.+b, 05.45.Tp, 89.75.-k, 87.19.La 
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Cortical networks exhibit diverse patters of activity, in- 
cluding oscillations, synchrony and waves. During neu- 
ronal activity, each neuron can receive inputs by thou- 
sands of other neurons and, when it reaches a threshold, 
redistributes this integrated activity back to the neuronal 
network. Recently it has been shown that another mode 
of activity is neuronal avalanches, with a dynamics sim- 
ilar to self-organized criticality (SOC) [J H H H, ob- 
served in organotypic cultures from coronal slices of rat 
cortex Je| where neuronal avalanches are stable for many 
hours |6(. The term SOC usually refers to a mechanism 
of slow energy accumulation and fast energy redistribu- 
tion driving the system toward a critical state, where the 
distribution of avalanche sizes is a power law obtained 
without fine tuning: no tunable parameter is present in 
the model. The simplicity of the mechanism at the basis 
of SOC has suggested that many physical and biological 
phenomena characterized by power laws in the size dis- 
tribution, represent natural realizations of the SOC idea. 
For instance, SOC has been proposed to model earth- 
quakes 0,11], the evolution of biological systems 0, solar 
flare occurrence , fluctuations in confined plasma [ill ] 
snow avalanches |l2j and rain fall [l3j ]. 

In order to monitor neural activities, different time 
series are usually analysed through power spectra and 
generally power-law decay is observed. A large number 
of time series analyses have been performed on medical 
data that are directly or indirectly related to brain activ- 
ity. Prominent examples arc EEG data which are used 
by neurologists to discern sleep phases, diagnose epilepsy 
and other seizure disorders as well as brain damage and 
disease 0,0. However, the interpretation of physiolog- 
ical mechanisms at the basis of EEG measurements is still 
controversial. Another example of a physiological func- 
tion which can be monitored by time series analysis is the 
human gait which is controlled by the brain [nj. For all 



these time series the power spectrum, i.e. the square of 
the amplitude of the Fourier transformation double loga- 
rithmically plotted against frequency, generally features 
a power law at least over one or two orders of magnitude 
with exponents between 1 and 0.7. Moreover, experi- 
mental results show that neurotransmitter secretion rate 
exhibits fluctuations with time power law behaviour |l7^ 
and power laws are observed in fluctuations of extended 
excitable systems driven by stochastic fluctuations (lif . 

Here we present a model based on SOC ideas and tak- 
ing into account synaptic plasticity in a neural network. 
With this model we analyse the time signal for electri- 
cal activity and compare the power spectra with EEG 
data. Plasticity is one of the most astonishing proper- 
ties of the brainjoccuring mostly during development 
and learning [UJ, |2fJ, |2l| , and can be defined as the abil- 
ity to modify the structural and functional properties 
of synapses. Modifications in the strength of synapses 
are thought to underly memory and learning. Among 
the postulated mechanisms of synaptic plasticity, the ac- 
tivity dependent Hebbian plasticity constitutes the most 
fully developed and influential model of how information 
is stored in neural circuits |U |U 0] ■ A large variety of 
models for brain activity has been proposed, based for in- 
stance on the convolution of oscillators [111 or stochastic 
waiting times [jfij] . They are essentially abstract repre- 
sentations on a mesoscopic scale, but none of them is 
based on the behaviour of a neural network itself. In 
order to get real insights on the relation between time 
series and the microscopic, i.e. cellular, interactions in- 
side a neural network, it is necessary to identify the basic 
ingredients of the brain activity possibly responsible for 
characteristic scale-free behaviour observed through the 
spectrum power law. 

In order to formulate a new model to study EEG sig- 
nals, we introduce within a SOC approach the three 



most important ingredients for neuronal activity, namely 
threshold firing, neuron refractory period and activity- 
dependent synaptic plasticity. We consider a simple 
square lattice of size LxLon which each site represents 
the cell body of a neuron, each bond a synapse. There- 
fore, on each site we have a potential vt and on each 
bond a conductance gij. Whenever at time t the value 
of the potential at a site i is above a certain threshold 
Vi > Dmax, approximately equal to — 5bmV for the real 
brain, the neuron fires, i.e. generates an "action poten- 
tial" , distributing charges to its connected neighbours in 
proportion to the current flowing through each bond 

v j (t + l)=v j (t)+v i (t) ' i f ) (1) 

where Vj(t) is the potential at time t of site j, nearest 
neighbor of site i, iij — gij{vi — Vj) and the sum is ex- 
tended to all nearest neighbors k of site i that are at 
a potential Vk < Uj. After firing a neuron is set to a 
zero resting potential. The conductances are initially all 
set equal to unity whereas the neuron potentials are uni- 
formly distributed random numbers between w max — 2 
and t> max ~ 1- The potential is fixed to zero at top and 
bottom whereas periodic boundaries are imposed in the 
other direction. 

The external stimulus is imposed at one input site in 
the centre of the lattice, and the electrical activity is 
monitored as function of time by measuring the total 
current flowing in the system. The firing rate of real 
neurons is limited by the refractory period, i.e. the brief 
period after the generation of an action potential during 
which a second action potential is difficult or impossible 
to elicit. The practical implication of refractory periods 
is that the action potential does not propagate back to- 
ward the initiation point and therefore is not allowed to 
reverberate between the cell body and the synapse. In 
our model, once a neuron fires, it remains quiescent for 
one time step and it is therefore unable to accept charge 
from firing neighbours. This ingredient indeed turns out 
to be crucial for a controlled functioning of our numerical 
model. In this way an avalanche of charges can propagate 
far from the input through the system. 

As soon as a site is at or above threshold w max at a 
given time t, it fires according to Eq. (1), then the con- 
ductance of all the bonds, connecting to active neurons 
and that have carried a current, is increased in the fol- 
lowing way 

9i j (t+l)=g ij (t) + Sg ij (t) (2) 

where Sgij(t) = fco£y(t), with a being a dimensionlcss 
parameter and k a unit constant bearing the dimen- 
sion of an inverse potential. After applying Eq. (2) the 
time variable of our simulation is increased by one unit. 
Eq. (2) describes the mechanism of increase of synap- 
tic strength, tuned by the parameter a. This parameter 
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FIG. 1: Total current flowing in one square lattice config- 
uration (L = 1000, a = 0.03, a t = 0.0001, tw = 6) as 
function of time in a sequence of several thousand stimuli. In 
the inset we show the asymptotic value of the percentage of 
active bonds as function of a for L — 100. The value of the 
parameters is at = 0.0001 and t) max = 6. 



then represents the ensemble of all possible physiological 
factors influencing synaptic plasticity, many of which are 
not yet fully understood. 

Once an avalanche of firings comes to an end, the con- 
ductance of all the bonds with non-zero conductance is 
reduced by the average conductance increase per bond, 
Ag = t ^9ij where Ni, is the number of bonds 

with non-zero conductance. The quantity Ag depends 
on a and on the response of the brain to a given stimu- 
lus. In this way our electrical network "memorizes" the 
most used paths of discharge by increasing their conduc- 
tance, whereas the less used synapses atrophy. Once the 
conductance of a bond is below an assigned small value 
at, we remove it, i.e. set it equal to zero, which corre- 
sponds to what is known as pruning. This remodelling 
of synapses mimicks the fine tuning of wiring that occurs 
during "critical periods" in the developing brain, when 
neuronal activity can modify the synaptic circuitry, once 
the basic patterns of brain wiring are established [2p| . 
These mechanisms correspond to a Hebbian form of ac- 
tivity dependent plasticity, where the conjunction of ac- 
tivity at the presynaptic and postsynaptic neuron mod- 
ulates the efficiency of the synapse 24]. To insure the 
stable functioning of neural circuits, both strengthening 
and weakening rule of Hebbian synapses are necessary 
to avoid instabilities due to positive feedback J2t|. How- 
ever, differently from the well known Long Term Poten- 
tiation (LTP) and Long Term Depression (LTD) mecha- 
nisms, in our model the modulation of synaptic strength 
does not depend on the frequency of synapse activation 

mum. 

The external driving mechanism to the system is im- 
posed by setting the potential of the input site to the 
value w m ax, corresponding to one stimulus. This external 
stimulus is needed to keep functioning the system and 
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FIG. 2: (Color online) Log-log plot of the distribution of 
avalanche size n(s) (L = 1000, a = 0.03 and 0.08, N p = 10, 
«max = 6) for the square lattice (lines) and the small world 
lattice (*, L = 1000, a = 0.05, N p = 1000, v ma * = 8) with 1% 
rewired bonds. The data are averaged over 10000 stimuli in 
10 different configurations. The dashed line has a slope 1.2. 
For a = 0.3 and random input site the slope is 1.5. 

therefore mimicks the living brain activity. We let the 
discharge evolve until no further firing occurs, then we 
apply the next stimulus. Fig.l shows the electrical sig- 
nal as function of time: the total current flowing in the 
system is recorded in time during a sequence of succes- 
sive avalanches. Data show that discharges of all sizes 
are present in the brain response, as in self-organized 
criticality where the avalanche size distribution scales 
as a power law 0, |3(| The strength of the parameter 
a, controlling both the increase and decrease of synap- 
tic strength, determines the plasticity dynamics in the 
network. For large values of a the system strengthens 
more intensively the synapses carrying current but also 
very rapidly prunes the less used connections, reaching 
after a short transient a plateau where it prunes very few 
bonds. On the contrary, for small values of a the sys- 
tem takes more time to initiate the pruning process and 
slowly reaches a plateau. The number of active (non- 
pruned) bonds asymptotically reaches its largest value 
at the value a = 0.03 (inset of Fig.l). This could be in- 
terpreted as an optimal value for the system with respect 
to plastic adaptation. 

Since each avalanche may trigger the activity of a high 
number of neurons, large currents flow through the sys- 
tem, therefore after N p stimuli the network is no longer 
a simple square lattice due to pruning, but exhibits a 
ladder-like pattern with few lateral connections. This 
complex structure constitutes the first approximation to 
a trained brain, on which measurements are performed. 
These consist of a new sequence of stimuli at the input 
site, by setting the voltage at threshold, during which 
we measure the number of firing neurons as function of 
time. This quantity corresponds to the total current flow- 
ing in a discharge measured by the electromagnetic signal 
of the EEC We have evaluated the size distribution of 
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FIG. 3: (Color online) Power spectra for experimental data 
and numerical data (L = 1000, a = 0.03, N p = 10, v maK = 6) 
for the square lattice (middle line) and the small world lattice 
(bottom line, L = 1000, a = 0.05, N p = 1000, u max = 8) with 
1% rewired bonds. Spectrum for the square lattice a = 0.3 
and leaky neurons = 0.01). The experimental data (top 
line) are from ref.[19] and frequency is in Hz. The numerical 
data are averaged over 10000 stimuli in 10 different network 
configurations. The dashed line has a slope 0.8. 

neural avalanches, that is the total number of neurons 
involved in the propagation of each stimulus. This dis- 
tribution exhibit power law behaviour, with an exponent 
equal to 1.2 ± 0.1, quite stable with respect to parame- 
ters (Fig.2). We have also simulated the brain dynamics 
on a square lattice with a small fraction of bonds, from 
to 10%, rewired to long range connections correspond- 
ing to a small world network [3lL I32I l33j . which more 
realistically reproduces the connections in the real brain. 
Fig.2 shows the size distribution scaling with an expo- 
nent 1.2 ± 0.1 for a system with 1% rewired bonds and 
a different set of parameters a, N p , t> ma x- Conversely, 
for the input site chosen at random in the system, the 
scaling exponent changes and becomes 1.5 ± 0.1 (Fig.2). 

In order to compare with medical data, we calculate 
the power spectrum of the resulting time series, i.e. the 
square of the amplitude of the Fourier transform as func- 
tion of frequency. The average power spectrum as func- 
tion of frequency is shown in a log-log-plot with the pa- 
rameters a — 0.03, N p — 10, o t = 0.0001, w max = 6 and 
a lattice of size L = 1000 (Fig. 3). We see that it ex- 
hibits a power law behaviour with the exponent 0.8 ±0.1 
over more than three orders of magnitude. This is pre- 
cisely the same value for the exponent found generically 
on medical EEG power spectra 0, |35| . We also show 
in Fig. 3 the magnetoelectroencephalography (similar to 
EEG) obtained from channel 17 in the left hemisphere 
of a male subject, as measured in ref.j3^, having the 
exponent 0.795. 

We have checked that the value of the exponent is 
stable against changes of the parameters a, w m ax, &t, 
and N p , and also for random initial bond conductances. 
Moreover, the scaling behaviour remains unchanged if 
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the input site is placed at random in the system at each 
stimulus. For a = the frequency range of validity of 
the power law decreases by more than an order of mag- 
nitude. Fig. 3 also shows the power spectrum for a small 
world network with 1% rewired bonds and a different set 
of parameters a, N p , w max : the spectrum has some devi- 
ations from the power law at small frequencies and tends 
to the same universal scaling behaviour at larger frequen- 
cies over two orders of magnitude. The same behaviour 
is found for a larger fraction of rewired bonds. 

In real systems neurons have a leakage, namely the 
potential decays exponentially in time with a relaxation 
time r, i.e. — — 7u(t), with 7 = 1/t. Leakage 

has been considered in our model and the same scal- 
ing behaviour recovered (Fig. 3). However for r < 10 
(i.e. for stronger leaking), the low frequency part of the 
spectrum appears to be frequency independent and the 
scaling regime is recovered at high frequencies with an 
exponent in agreement with previous results. 

In the mature living brain synapses can be excitatory 
or inhibitory, namely they set the potential of the post- 
synaptic membrane to a level closer or farther, respec- 
tively, to the firing threshold. We have introduced in our 
model this ingredient: each synapse is inhibitory with 
probability pi n and excitatory with probability 1 — Pi n . 
We have studied the power spectrum for a range of value 
of pi n . For a density up to 10 % of inhibitory synapses 
the same power law behaviour is recovered within er- 
ror bars. For increasing density the scaling behaviour is 
progressively lost and the spectrum develops a complex 
multi-peak structure for pi„ = 0.5. Furthermore, the 
size distribution exhibits an exponential behaviour even 
for very small densities of inhibitory synapses. These 
results suggest that the balance between excitatory and 
inhibitory synapses has a crucial role on the overall be- 
haviour of the network, similarly to what can occur in 
some severe neurological and psychiatric disorders |36|. 

The stability of the spectrum exponent suggests that 
an universal scaling characterizes a large class of brain 
models and physiological signal spectra for brain con- 
trolled activities. Medical studies of EEG focus on sub- 
tic details of a power spectrum (e.g. shift in peaks) 
to discern between various pathologies. These detailed 
structures however live on a background power law spec- 
trum that shows universally an exponent of about 0.8, 
as measured for instance in refs. |34| and |35| . A sim- 
ilar exponent was also detected in the spectral analysis 
of the stride-to-stride fluctuations in the normal human 
gait which can directly be related to neurological activ- 
ity 0| . Our simple model is based on SOC ideas: the 
threshold dynamics ensures time scale separation (slow 
external drive and fast internal relaxation) . This dynam- 
ics leads to criticality and therefore power law behaviour 
Q. However the new ingredients of the model, namely 
the plasticity of the synapses may be at the origin of the 
new observed exponent. This work may open new per- 



spectives to study pathological features of EEG spectra 
by including further realistic details into the neuron and 
synapsis behaviour. 
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